{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__title__      = ML in Action Chapter8 Code（预测鲍鱼年龄）\n",
    "__author__     = wgj\n",
    "__createDate__ = 2018-11-03 21:20:19"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\"\"\" 加载数据集 \"\"\"\n",
    "def loadDataSet(fileName): \n",
    "    numFeat = len(open(fileName).readline().split('\\t')) - 1 #  数据集列数（末列为真值）\n",
    "    dataMat = []; labelMat = []\n",
    "    fr = open(fileName)\n",
    "    for line in fr.readlines():\n",
    "        lineArr =[]                                                 # 保存一行训练数据\n",
    "        curLine = line.strip().split('\\t')                # 移除每行首尾空格并用tab分割数据\n",
    "        for i in range(numFeat):                        # 遍历所有列\n",
    "            lineArr.append(float(curLine[i]))\n",
    "        dataMat.append(lineArr)\n",
    "        labelMat.append(float(curLine[-1]))\n",
    "    return dataMat,labelMat\n",
    "\n",
    "\"\"\" 标准线性回归 \"\"\"\n",
    "def standRegres(xArr,yArr):\n",
    "    \"\"\" \n",
    "    输入：1、xArr——输入变量\n",
    "                2、yArr——目标值\n",
    "    输出：ws——回归系数\n",
    "    \"\"\"\n",
    "    xMat = mat(xArr); yMat = mat(yArr).T\n",
    "    xTx = xMat.T*xMat\n",
    "    if linalg.det(xTx) == 0.0:\n",
    "        print(\"This matrix is singular, cannot do inverse\")\n",
    "        return\n",
    "    ws = xTx.I * (xMat.T*yMat)\n",
    "    return ws\n",
    "\n",
    "\"\"\" 局部加权线性回归 \"\"\"\n",
    "def lwlr(testPoint,xArr,yArr,k=1.0):\n",
    "    \"\"\" \n",
    "    输入：1、testPoint——测试点\n",
    "                2、xArr——输入变量\n",
    "                3、yArr——目标值\n",
    "                4、k——权值\n",
    "    输出：testPoint * ws——预测结果\n",
    "    \"\"\"\n",
    "    xMat = mat(xArr); yMat = mat(yArr).T\n",
    "    m = shape(xMat)[0]\n",
    "    weights = mat(eye((m)))                  # 创建对角矩阵\n",
    "    for j in range(m):                             # 计算高斯核对应的权重\n",
    "        diffMat = testPoint - xMat[j,:] \n",
    "        weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2))\n",
    "    xTx = xMat.T * (weights * xMat)\n",
    "    if linalg.det(xTx) == 0.0:\n",
    "        print(\"This matrix is singular, cannot do inverse\")\n",
    "        return\n",
    "    ws = xTx.I * (xMat.T * (weights * yMat))\n",
    "    return testPoint * ws\n",
    "\n",
    "\"\"\" 逐点调用lwlr() \"\"\"\n",
    "def lwlrTest(testArr,xArr,yArr,k=1.0): \n",
    "    \"\"\" \n",
    "    输入：1、testArr——测试点集合\n",
    "                2、xArr——输入变量\n",
    "                3、yArr——目标值\n",
    "                4、k——权值\n",
    "    输出：yHat——预测结果\n",
    "    \"\"\"\n",
    "    m = shape(testArr)[0]\n",
    "    yHat = zeros(m)\n",
    "    for i in range(m):\n",
    "        yHat[i] = lwlr(testArr[i],xArr,yArr,k)\n",
    "    return yHat\n",
    "\n",
    "\"\"\" 计算回归误差（平方和）\"\"\"\n",
    "def rssError(yArr,yHatArr): \n",
    "    \"\"\" \n",
    "    输入：1、yArr——目标值\n",
    "                2、yHatArr——预测结果\n",
    "                （计算元素平方，输入需均为arrays类型）\n",
    "    输出：回归误差\n",
    "    \"\"\"\n",
    "    return ((yArr-yHatArr)**2).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "测试局部加权线性回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4VGX2wPHvSaPXJCBSjFKkKaBRCMqKihJXQUVERCwrgmXXAlbQJQbEtWP5sSKKa1sLKLKIFBEFlAzBoPReJRQJgQQChITk/P6YiYQwSSYwyZScz/Pkyczcd+49l4Qzb8593/eKqmKMMSa4hPg6AGOMMd5nyd0YY4KQJXdjjAlCltyNMSYIWXI3xpggZMndGGOCkCV3Y4wJQpbcjTEmCFlyN8aYIBTmqwNHRUVpTEyMrw5vjDEBacmSJXtVNbq0dj5L7jExMaSkpPjq8MYYE5BEZJsn7awsY4wxQciSuzHGBCFL7sYYE4QsuRtjTBCy5G6MMUHIkrsxxgQhS+7GGBOELLkbY0wQsuRujDEFHA6Ij3d+D3A+m6FqjDF+JzERZs92Pp41y7exnCZL7sYYUyAh4cTvAcySuzHGFIiLC/geewGruRtjjLf4Uc3ekrsxpnJyl4hPNzkX1OwTE70T42mwsowxpnJyJeIDB4+w/oxz2H/bXWR8+TWN1u6mS+IowmbNLPs+/ahmbz13Y0zwKKnnXWjbsbx8/jfoKQbd9xaxlwyjb8u+DF6cxePNrmJg/zF073APn9z2GNk/J5XtWHFxzsSemOjdvwhOgfXcjTHBo6ShjImJ5M75nq/rt2Fcp2y2pR/mzGZtuD1KufS7yUTdPZA6L45h9a4DjO96C880vZw3vk5lcP4mBnQ+i5pVXOnS4YChQ2HNGjhw4ORjuYvBB0MsPU7uIhIKpAA7VPU6N9v7Ac8CCixT1QHeCtIYYzxSUA7p08fZU05IgLg4cvPy+fqup/i/Nnfwe5U6tKsSxoSLa9Jjwr8ISUiAwe8435eYQbP1yfQMy8TRpgv/juvH8zPWMu7HTdzVNYa7usZQLzERkpOd7SMjTy7BuCvN+KBcU5ae+8PAGqB20Q0i0hIYDlyiqvtFpIGX4jPGGM8VDGWMj4fZs8mREKYkjmfcvI1s33eE885pxntXtuTKNg2Qa645uTc9diwMHYoAXR+9h65xcSzdnsG/f9zIG3M38O5Pm7ntpscZfPgYDXKynO3j4tzH4GMeJXcRaQJcC4wBhrlpMhgYp6r7AVR1j9ciNMaYMsr550i+jGzHuFY92DFlBR3qhpK4+Rsu730H0rahs5G73nRcHNSte3zEy6xZdGxalwl3xLJu9k+8/dViJua05sNuj3JzbBPubdWcZp4E5MdlmdeBJ4BaxWxvBSAiC4FQ4FlVPekMRGQIMASgWTOP/kmMMaZM9h/K4dZfjrG26RV0iKrNcz1a0v3O3khyMvy+FBYtcjYsroddTAnl3LFjeH32bIb1upnxf/snk1NS+fyX7fRo04BGdapRJTyEEBFCRWjeoAaXtIiiQa2qJe6zPJWa3EXkOmCPqi4Rke4l7Kcl0B1oAvwkIu1VNaNwI1WdAEwAiI2N1dOI2xhj3Hr83Xls3nmE8Z3r0rNPV0SkbDsoJek3Gz6U5+PO4+ErW/LeT5v5dvkukjalc/RYPqrKsXxFXdmt9Rm1uL5jYwZcHEudCi7VeNJzvwToLSJ/BaoCtUXkE1UdWKhNKrBIVXOBLSKyDmey/8XrERtjTDF+3rCX73fn8uSCj4lfngU3XebcMHasszRyOj3nIkm/Ye2qPF0/k6fnv/znhVuA/Hxl1c4D/LQxjXlr03hx1lrenLuBfrFNGHbVudSpHn46p+gxUfW8A+3quT9WdLSMiMQDt6rqnSISBfwGdFTV9OL2FRsbqykpKacWtTHGuHHXfxazduNu5r87hCrPjYIhQ8r3gK4Lt/TsWWwtffXOA7y/cAtTf9tBo7pV+feACzmvSZ1TPqSILFHV2NLanfIkJhEZJSK9XU9nA+kishr4EXi8pMRujDHetivzCAvWp9FvxxKqpP0BU6aU38EKJiX16eNM7O7+InA4oF072rZtxiv7k5l0Xxx5ecrzM9ZQlk71qSrTJCZVnQfMcz0eWeh1xTmKxt1IGmOMKXdTft1BvkLfgVfB7uTyvXhZMPolI8M5uqa4NqtXOx+PGMEFQ4YwPa4qua+8jJz32MlDKL3MZqgaYwKeqjIpZTtx50TSrEcX6FHOFy8LPjgyMoof4piQANu3w44d8PzzANR/YbSz/bFD5T4k0taWMcYEvMVb9rEt/TD9LmpSMQcsuLg6dqyzLFMwI7boOjOrVjk/AIYMcW7LyIDOnStkSKT13I0xAW9SSiq1qoQR365RxR64yIxYoPgeecGyBT17lntJBiy5G2MC3MHsXGas2MUNnRpTLSLUN0F4MkmpgicyWXI3xgS0b5fv4khuHv1iK6gk444n68lU8JozVnM3xgS0SSnbadmgJh2bFjNqpZKy5G6MCVgb9xzk198z6BfbtOzLDAQ5S+7GmIA1OSWVsBDhhk6NfR2K37HkbowJSDnH8vly0Rau2LeR6JW/+jocv2MXVI0xAWn2qt2k5ygDvvsQts7wixtk+BPruRtjAotrXZdP56ygSfUQ/tIqukLXSQ8UltyNMYElMZHNi1fg2HuMW7u1JGTWrAqZFBRorCxjjAksCQlMnDifiBC42Zdj2/2cJXdjTED5o10nJjfMoG9sk+O3sTMnsbKMMSag/PvHjeSpcv9lzX0dil/zOLmLSKiI/CYi00to01dEVERKvUuIMcaU1S9b9/HRom0MuLgZTetX93U4fq0sPfeHgTXFbRSRWsBDQPLpBmWMMUUdzjnG45OX0aReNZ6KPHDyErvmBB4ldxFpAlwLvFdCs9HAS0C2F+IyxpgTvDRrHVvTD/Ny3w7UGDPKucRuYqKvw/JbnvbcXweeAPLdbRSRTkBTVS22ZONqN0REUkQkJS0trWyRGmMqLcemdD5I2spdXWPock6kc1x7cfcuNYAHyV1ErgP2qOqSYraHAGOBR0vbl6pOUNVYVY2Njo4uc7DGmMpnz4Fshk1aSkxkdZ6Mb+18sWD5XBvfXixPeu6XAL1FZCvwOXCFiHxSaHstoD0wz9WmCzDNLqoaY07XkZw87vkohcwjuYy77QLf3YwjAJWa3FV1uKo2UdUYoD/wg6oOLLQ9U1WjVDXG1WYR0FtVU8oraGNM8MvPV4ZNWsqKHZm82aka7e6+xS6glsEpj3MXkVEi0tubwRhjKjnXujE4HLz83TpmrtzN039tQ493/mUXUMuoTDNUVXUeMM/1eGQxbbqfblDGmEoqMRFmz2ZS/Ta83awHAzo3Y9ClZ1f4/UeDgS0/YIzxHwkJOGo2ZsRZPejWIorE3u2cd1iq4PuPBgNL7sYYv7G5xXncd94txNSqwv8NuIDwUFsh5VTZv5wxxi/k5ytDv1hKaIjw/p0XUadauK9DCmiW3I0xfuGb5TtZlprJiL+2oVmkrRtzuiy5G1NZFRqZ4mtHj+Xx8ux1tG1UmxvtZtdeYTV3Yyor18gUwOcXKz9M2krq/iN8Muh8QkPEp7EEC+u5G1NZlbQ+SwX26tfsOsCr363nytYNuLRlVLkfr7Kw5G5MZVOQuKH49VkKevXlPGko6+gx/v7fX6lTLZwX+55frseqbCy5G1PZFCTuXr2K75kX7dW768mfZu9eVRkxZQVb0w/x5vkRRPW93i/q/8HCau7GVDJZI0by264jbA6rTec77qf1R2+f3HsvOmnIXX3+NGv2n/+ynWnLdvLY1a3o8tKDflP/DxaW3I0JRg6HM/kmJPyZuLNz8/j3vE2Mn59JzjVP/Nn04v+uoF9EU9o0qsU5UTXdr7zobvq/J0sCuIkDYPXOAyRMW0W3llE80L0FVLPlBbxNVNUnB46NjdWUFFs40phyER/v7An37AmzZvHj2j0kTFvF7/sO06vDmfSLbcJZ29Yx+8Nv+ajlX9h++Ph9eM6sU5Wzo2twTlRNWjeqRY82DWlYu6pzo8NB5nP/4oc7hzE/vw6NDu3jpmnv0eKph9zX7ovEAc46e++3fibr6DFmPNyNqJpVKuJfJGiIyBJVLXVJdUvuxgQjV485Y/g/Gb70EDN35dK8Zgij+19E1xYnjkjJy1fW7T7Ilr2H2JyWxea9h5xfaVkczD4GwIVn1SPunEg2fD2bH2o1Izc0nPo1IjiQdYR8hZv3r+WRVx6kUZ1qbuMo6Llv33eYxyYv45et+/h0cBfnXZVMmXia3K0sY0wwWrGCzeu3c/ucvezJER7/6RMG180i4pkZx9u4Em9oQgJt4+Joe2btE3ahqmxKy2Lmit3MXLmb//txI43OaMnA7b9yXZ9udLqmG/sXJDHuox/5uGEnpr48j79dcjb3d29+fOkAV+0+L1/5z0+befW79YTk5/HS1jl0+aM+nGN3Uiov1nM3JghtbH4e/a95Eg0NYeLqL+l4cCeMHXti6cRNyeQERXrdeflKaPIitzX07fsO89qc9UxduoPaVcP5xzlh3P7pqxx6cgSO+mczYcFmlqdmcvm50Tw35SUaT/+q+OOaEnm95y4ioUAKsENVryuybRhwD3AMSAPuVtVtZQvZGOMN2bl5DBkwGg4c5otFE2ixeL4zkRatiRdcvOzTx5noiyTsoqNhQpMXOYdPpqf/+VqBpvWrM/aWjgzudg4vzV7LmFVpvNxuMDmzMoDfaFy3Gm/070jvDmcirR+F3Cy7eFreVNWjL2AY8Ckw3c22y4Hqrsf3A1+Utr8LL7xQjTHeN3bOOj3ryem6YP0e1aQk1Z49nd+L07OnKji/F1b4vUlJqpGRznaRkSXvT1WTps3XxL+N1nc/mquLNu3VYwsXlh6H8QiQop7kbI8aQRNgLnCFu+RepG0nYGFp+7TkboyXJSXp1l43a8vh0/Ufn/5apvd5/AHgQWIv8f1FP0BMmXma3D0ty7wOPAHU8qDtIGCmh/s1xniJJiaSUCuOiNwcnrm2jedv9OQuR4XHtLsb8lgau01ehSt1+QERuQ7Yo6pLPGg7EIgFXi5m+xARSRGRlLS0tDIHa4wp3ndDhjOveSxDGxymYb8bvDuVv+AD4FQSuzfeb8rMk7VlLgF6i8hW4HPgChH5pGgjEekBPA30VtWj7nakqhNUNVZVY6Ojo08jbGNMYfn5ymvbhObRNbhz7icVsuhXqfxovfjKqNSyjKoOB4YDiEh34DFVHVi4jYh0At4B4lV1TznEaYwpwXdvfMy6PyJ5o8E+whJGAur7EogfrRdfGZ3yJCYRGYWzsD8NZxmmJjBZRAB+V9Xe3gnRGFMSVeWtVQc5OySb6yY+DWlp/pFMrc7uU2VK7qo6D5jnejyy0Os9vBqVMcZjP6zdw6qoGF6e/y6hY8b4OpzjPLlQa8qNLT9gTABTVd78YSNN6lXjhp+nQKjdosE42W+CMQHspw17WbY9gwe6tyDcErspxH4bjAlQqsqbczfQqE5Vbrqwsa/DMX7GkrsxAcqxOZ2Ubfu5v3tzqoS5ucGGqdQsuRsToN6au5EGtarQL2+njSc3J7ELqsYEoJSt+3BsTueZa9tQ9bkHbDy5OYkld2MC0Js/bCSyRgS3dT7LxpMbtyy5GxNgls74iQXrD/Bk22rOm1nbeHLjhtXcjQkwr01fQb3Dmdz++Vhfh2L8mCV3YwLI4i37WFD7LO7fv4Ka/xzh63CMH7OyjDEBQlV55bt1RNeqwu0Tn4MIG/5oimc9d2MCxMKN6Szeso9/XN7CWWs3pgSW3I0JAAW99sZ1q9H/4qa+DscEAEvuxgSAH9buYen2DB66soXNRjUeseRujJ/Lz1de/W49MZHV6XNBE1+HYwKEx8ldREJF5DcRme5mWxUR+UJENopIsojEeDNIYyqzWat2s3rXAR7p0cpWfjQeK8tvysPAmmK2DQL2q2oLYCzw4ukGZkyl53CQF38Nr/1vKS0b1KRXhzN9HZEJIB4ldxFpAlwLvFdMk+uBD12PvwSuFNf99owxpygxkSk7ctmYlc+wq1oRGmL/pYznPB3n/jrwBFCrmO2Nge0AqnpMRDKBSGDvaUdoTCV1cMRIXvxfKp3qhdKz3Rm+DscEmFJ77iJyHbBHVZeU1MzNa+pmX0NEJEVEUtLS0soQpjFe5nD4/TK5b647QnpYdZ79/h1Ckhf5OhwTYDzpuV8C9BaRvwJVgdoi8omqDizUJhVoCqSKSBhQB9hXdEeqOgGYABAbG3tS8jemwiQm+vUyuat3HuD9jYe5ZfkcOnw/FUKP+GWcxn+VmtxVdTgwHEBEugOPFUnsANOAOwEH0Bf4QVUteRv/5cfL5OblKyO+XkHdKqE8dXA5dO7sl3Ea/3bK46pEZJSI9HY9nQhEishGYBjwlDeCM6bcFCyTGxfnfruvyjYOB5/e/jhLt2fwTPvq1K1ZFcaOLT5OY4pRpoXDVHUeMM/1eGSh17OBm70ZmDE+5YuyjcPBnn4Deemm57nk4HZueH+qX5eOjH+zGRGm8vGkV56QAD17nlgOKe/efGIioy7oy9GwCEbf0B5xF4MxHrIlf03lU7RX7nA4X0tIOF7+cHd3o3Luzf9433CmL8piaOtqnHPVpeV2HFM5WHI3lc7RZ0Yyp24LIvrdTMPtGZwx5hWiZ39HCJA3YyZHcvM4mptH/RoRnDAXrxwvwu7IOMKjvx6i5ZF07mvRwuv7N5WPJXdTqRzIzuWWJXmsibkGFmfB4oXQ/m6qtr2d/LBwckbM+LNt47rV6NXhTHp3OJM2jWq5nczhDVlHj3H/J0vIPZLN+E9GUGVde+uxm9Nmyd1UGvn5ytDPl7Lhj4OMG3ABTetXY3dmNrsys9mafogqYaFUCw+lWkQIISL8tGEv7/60mfHzN9E8uga9l8+l/8JfaJiY6JXkq0lJfPvmp4xu/Vf2HFXeiYuk+Yr2VmM3XiG+Go4eGxurKSkpPjm2qZxem7OeN+duYNT17bgjLub4Bnc1d5d9h3KYsWIX3yzbyeIt+2iQk8UX1zQmpqAmXho3+85akMSM8ZP5os65LKnTlHaH9zD6sRu4oFk9L52pCWYiskRVY0trZ6NlTPBzOPiu3/28OXcDfS9swu1dzjpx+9ChzgulQ4eeNCKmfo0IBnY5iy/ujWPGw93IqVuPvy07xv5DOZ4du+AibGIi2bl5vD1vE3Hf7OaJZj3YF16d0dvnMa3/uZbYjddZcjdBL/X513i08eWcn7aZ5xodosQFSwsl46LaNKrNu3fEsmP/Ee79eAlHj+WVPjwyIYH8nj2ZevdTXPnqfF6ctZbOjWvy1frJ/NC/Bbf/92VCu3b10pkac5zV3E1Qy89XnrjsHvJ3HmTcV2OouqPjyfXysWOPl04KFFP3jo2pzyv9OvDQZ7/x5JfLGftBIlLC8Mi1Z7fj8WufYcWvmbRvXJuXbz6frs2j4JGrvXWKxrhlyd0EtVfnrCNp7zH+dVEkTVd3dJ+0i45pL+Viae/D29i+08HLxNHslscZBm73+8vWfdz9wS9UDQ9l7C0duL5DY0JsTXZTQawsY4LWl0tSGffjJvpf1JT+/buXvJaMpxwO6NWLBz4eQ7/0Vby5Lpv3n3n7pP3OXfMHA99LJrpWFab+/RJu7NSk+MQeAMsPm8Bjyd0EF1eiTP5mAcOnLKdr80jnVH5v3RgsMRHS05HISMbcdSk92zVk1PTVvDl3AwUjz75aksqQj5dw7hm1mHxhOI3731hy4i6hzm/MqbKyjAkuiYlsTV7Ove330LRRFG/fdqF3bypdaJZqeFwc47rk88RXy3ltznq2fDENadaMKVnVuSQqjHcGd6Hm9deVvmSBHy8/bAKXjXM3wcE1njzz+pu4cV019tWqz9SHLyMmqka5Hzo/X3lx0GgmRnUgPD+PW5fO5MmIHVSZOaPEMfTGnAob524ql8REdPZsHlmew/Ya9ZlwdxdnYq+AenZIiDB8yFWsWvkOv7Xcy8jw7VQZ+U/nxtLWjTemnFhZxgSHhAQmRnfkxzoxJF7blovPru98vaLWZY+LcyZ066UbP+HJDbKrishiEVkmIqtE5KSrPiLSTER+FJHfRGS5636rxpQPN73x5U3b8GKzv3B124bcEVdoBmpFroluF0aNH/Gk534UuEJVs0QkHPhZRGaqauHbsT8DTFLVt0WkLTADiPF+uMZwUm/8YHYuD372G9E1q/BS3/NPHBnjbl328uBwQEaG3e/U+A1PbpCtQJbrabjrq+hVWAVqux7XAXZ6K0BjTlJodImq8vTXK0ndf4QvhnShbvUI38SUmAjJyc6/EqwkY/yARzV3EQkFlgAtgHGqmlykybPAdyLyIFAD6FHMfoYAQwCaNWt2iiGbSq9Qb3xyynamLdvJo1e1Ijamvu9isuGMxs+UaSikiNQFvgYeVNWVhV4f5trXqyISB0wE2qtqfnH7sqGQ5nRt/O5nes3ZS8eG1fhk2NWE2tR+UwmUy1BIVc0A5gHxRTYNAia52jiAqkBUWfZtTFlk5+bxj283Ue3oYV6fN8ESuzFFeDJaJtrVY0dEquEsuawt0ux34EpXmzY4k3uad0M15rgx365hbbUoXt01n4ZPP+brcIzxO57U3BsBH7rq7iE4R8VMF5FRQIqqTgMeBd4VkaE4L67epb6a+mqC3qyVu/l40TYGdzuby194x9fhGOOXPBktsxzo5Ob1kYUerwYu8W5oxpzsYHYuz0xdSfvGtXm8Z2vnizbF35iT2AxVE1De+HgB6QezeX/+OCJiw5zJvKJmoRoTQCy5m4Cx/o+D/GfjYfovm83530+F0CPOZG7DEI05iSV3ExBUlWenraJmKDy+dCq0bXs8mVfULFRjAoitCmkCwowVu0nalM5j2+ZT/48d0LSp1deNKYH13I3fO3T0GM99u5q2jWozoHcfSF9pJRhjSmHJ3fi9cT9uZFdmNm/d2onQmPpWgjHGA1aWMX5tc1oW787fRJ99a4jdtc7X4RgTMKznbvyWqpL4zWqqHMvhqU/HwO+drddujIes52781vdfL2D++jQeiT5Mg0ttnXRjysJ67sYvHT2Wx+ifdtDyYAZ3Lp1uPXZjysh67sYvfezYxu9V6vBM5lLCE/7p63CMCTjWczd+Z9+hHN6Yu4HLWkVz2QsTfB2OMQHJkrvxLw4Hb773I4eiO/D0tW18HY0xAcvKMsavbHrhDT6p345b966gVcNavg7HmIBlyd04l8yNj3d+97F/XX43Vcln6MBuvg7FmIDmyZ2YqorIYhFZJiKrRCSxmHb9RGS1q82n3g/VeE3RZF6wZG6i2x9thUmatoDvd+fyQPs6RHW32wMYczo8qbkfBa5Q1SwRCQd+FpGZqrqooIGItASGA5eo6n4RaVBO8RpvKLr+uR8smZuXrzz3/WYa5+Zz9+dfwZ1X+CwWY4KBJ3diUiDL9TTc9VX0FnqDgXGqut/1nj3eDNJ4WdFk7gdL5k75NZXV1aN5Y+ssqo58xqexGBMMPKq5i0ioiCwF9gBzVDW5SJNWQCsRWSgii0Qk3tuBGi+Ki3Mm9sREv6izH845xsuz19GxaV16f/amLeVrjBd4lNxVNU9VOwJNgItFpH2RJmFAS6A7cCvwnojULbofERkiIikikpKWlnZ6kZvT42mdvQIutk5YsJk9B4/yz+vaICLldhxjKpMyjZZR1QxgHlC0Z54K/E9Vc1V1C7AOZ7Iv+v4JqhqrqrHR0dGnGLIpE3fJ2eGAjAzo3BkdOZID2bnk5uW7b1/OF1t3Z2bzzvzNXHtmBBfeO8Av/pIwJhiUWnMXkWggV1UzRKQa0AN4sUizqTh77B+ISBTOMs1mbwdrToGbm0fvfv4VPqjamtkde7BjRiY5074jIiyENmfUov3Sn+m25QA9ExORCrjY+sLMNeTlK099/67d5NoYL/JktEwj4EMRCcXZ05+kqtNFZBSQoqrTgNnA1SKyGsgDHlfV9HKL2niuSHL+Ye0fDOs4iIM5eXRrWIWr2zQlskYEe7NyWJ6awbQzzue/N55HbP0wnt2RSfuyXGx1OJwfJgkJzrp50edF2n7x1mSmNruSh3cn0/S6HpBzwFZ+NMZLxDkYpuLFxsZqSkqKT45dGeUuTOKliXN5t8EFtGlUm3EDOnFOdM2T2uXnK1/+mspLs9aSnpVD/70refqObtT8S1f3Oy6cwAv+SujZ0/mBEB9/4vNCVt54O32a9+HiXev48LOnCb34Ili0yP0xjDF/EpElqhpbWjuboVoJHMnJ47bPVvJugwsYuHc5Xz/QlXM2rnB7oTQkeRH9nhnMD3+pzqA9vzGpfhv6frWB1P2H3e+8cE0+IcGZyAv/tVD4Oc4bcMxetZvbz7uVqPyjvLF6CqGaX16nbkzlpao++brwwgvVlL/cY3l6938Wa8yT03VK/4dUk5KcG3r2VAXn98IKv56UpAtuukfbP/2tXjh6ji7ekn7yAZKS/mxbmq17s/Su95P1rCen69WvzdctaVller8xRhVnObzUHGtlmSCmSUk89Z+f+SKyHc/d0J6BXc46vrG4erib1zfuOcigD1NI3X+EoT1acn/3FoSGeDZkUVVJ3X+EyUtSGT9/ExGhITzSoyV3do0hPNT+cDSmrDwty1hyD1YOB5899grDu93N33f/wuMfPHtauzuYncuIr1fyzbKddG0eydhbOtKwdtWT2mUeyWXZ9gyWur6Wbc8g/VAOANd3PJMRf21Dw1W/FX+h1RhTIkvuldyKG2/npuZ96Lx7PR8+0I2QrsVcEC0DVWVySiojp62kRkQYr5wXweXj/4WOHMncus1564cNLEvNBEAEmtcIoeO2lXS8PJaLr4w9voRvCRdajTEl8zS52806glB61lHu7TiA6MwMXr/nUq8kdhwOJDGRfgkJXPDgpfzj09/426KD9Alrx6bP1rCsxn5iIqvz2NWt6Ni0Huc3rUPt7t0gORlWdIYBhUbC+MFCZcYEOyt6BpncvHz+/umvpOfCO1c1IfKF0d6Z9VloVEyLBrWY+vdLuD2mClPbdSeDMF5c+w1zulbhH1e05NKWUdSuGn78vQcPnjgyp2DsvJVkjCk31nMPMv+5zP9OAAAPsklEQVSasZZFm/cx9pYOtL+/r7PnnJHhfgx5SZOMiirS264aHsro+3owsv15hK1aiQBkbz2xzDJ2rHP/GRk2+9SYCmbJPVg4HEx583PeP+tq7r7kbG7s1KT097hZmqBYxcxUDU/d7nwQGnpymaXgPYU/RIwxFcKSe6BzJc6Vx6oyvMMddDmYyvD69ZxlkLvvhrp1i0+q3qh9v/QSjBgBzz9ffO/fD9aLN6aysdEygS4+nv3zk7ju3rchJIRpvZs56+w2GsWYoGSjZSoJHTmSJyc62BNRiymX1yOy+yVQxUajGFPZ2WiZAPd5WBO+i27NE/M+4Lw3xjhfrMjRKBVwMw9jTNlZzz2A/XEgm9HTV3NpdBiD6h/2TU+9LBdljTEVxpJ7AHtl9jpyc/N4bv5HhPhqKr9NSDLGL3lyJ6aqwAKgiqv9l6rq9n+yiPQFJgMXqapdLS1HK3dk8uWvqdzzx6/EfDMJcjJ903O2kTDG+CVPeu5HgStUNUtEwoGfRWSmqp4wK0ZEagEPAcnlEKcpRFUZPX019apH8I+7roCdC63nbIw5QakXVF1LCGe5noa7vtyNnxwNvARkey88484Pa/eQvGUfQ9d+R53wEJvKb4w5iUejZUQkVESWAnuAOaqaXGR7J6Cpqk4vhxhNIfn5ymtz1tPsaCb9P3dN7zfGmCI8Su6qmqeqHYEmwMUi0r5gm4iEAGOBR0vbj4gMEZEUEUlJS0s71ZgrtdmrdrNq5wEe7nIm4Vf1sHKMMcatMo1zV9UMYB4QX+jlWkB7YJ6IbAW6ANNE5KQZVKo6QVVjVTU2Ojr6lIOurPLylbHTltE8ex83NI2wcowxpliejJaJBnJVNUNEqgE9gBcLtqtqJhBVqP084DEbLeN905fvZP3BPN6a/S6ha6vbKBVjTLE8GS3TCPhQREJx9vQnqep0ERmF80at08o1QgPAsYVJvD55M61r1+Tas6pbOcYYU6JSk7uqLgc6uXl9ZDHtu59+WKaoKeMms6VZD95Z9y0h1mM3xpTC1pYJADnH8nmz9dWcd/gPrn54oK/DMcYEAFt+IABMStlO6uF8Rv/9OuTcBr4OxxgTACy5+7nsn5P4vynbufCMunRvZSOMjDGesbKMn/ts/FR2R9Tk0cWTERFfh2OMCRCW3P3YkZw8xjXvTpeDqXR9dJCvwzHGBBAry/ixD5K2sveo8vbjN0NMfV+HY4wJIJbc/VTmgoWMn76byxvX4iJL7MaYMrKyjJ+a8P4cMsOq8tiCj3wdijEmAFly90Pb0g/x7pkXcf2+dbR74u++DscYE4AsufuhUR8tJDw3hxEdazuX9LWbTxtjyshq7n5m7po/mPtHLiMWfEzD/yyA9HTnBltywBhTBtZz9yPZuXkkfrOaFjVD+Fu9bHj+eejZ0xYJM8aUmfXc/cjb8zbx+77DfHpPZ8KfmeF8ccgQ3wZljAlI1nP3Bw4Hy2+8g3E/bOD6jmfStUVU6e8xxpgSVM7k7nBAfLzfXKg8MmoMj0R1JfpIJqM+ftZv4jLGBK7KmdwTE2H27JNvLu2jpP/89Y+wObIpr66dRp0Z0+ym18aY01ZqcheRqiKyWESWicgqETkp84jIMBFZLSLLRWSuiJxVPuF6SUICdO4MGRknJvLikn5xvPBh8M2ynXy89Sj3XHo2XRMetguoxhiv8OSC6lHgClXNEpFw4GcRmamqiwq1+Q2IVdXDInI/8BJwSznE6x1xcaQdVV5vEMva8Q5qOg5S84wo6tzwGAOqRdP+iQeKf6/D4Uz+CQnHPwzglIYqJm9O5/Evl3FRTD2eiG8NYSE25NEY4xWl9tzVKcv1NNz1pUXa/Kiqh11PFwFNvBqll81auZur/zKUSedfRVhuDhlbUlmz6wDTduVxY9tbGZ/TgPx8df/mwr37hATPe9pFevk/bUjjzomLaJyxh7ffeYSIS7tard0Y4zUeDYV03Rx7CdACGKeqySU0HwTMLGY/Q4AhAM2aNStbpKfL4SBr9PMkXvcQk3/PoX2dKrz+0Qha1BB47z2IiyPjcA7Dp6zghZlrmb8ujddu6UCjOtVO3E9BIk9IgLg4z3raDgf06vXnhKQpz7/HU1NWcE7WXj55/2GiDmc62yUmWs/dGOMVolpMD9VdY5G6wNfAg6q60s32gcA/gMtU9WhJ+4qNjdWUlJQyhnvqlt94B39veBk76jTg/nqHePi1h4l4bvRJ48hVlclLUnl22irCQ0N4oc95XHNeo7If0E355kiDRjyX+BH/3XqUzmfXZ3zrfOo9+SgcPAi1asHYsc4PDGOMKYaILFHV2NLalWm0jKpmAPOAeDcH7AE8DfQuLbFXGFcp5H+T53Fzm1vIr16dSd3q8PgL9xGRtgeeeOKkt4gI/WKb8u1D3YiJrM79//2VJ75cxqGjx8p2AbVI+WbVDbfR68GJ/HfrUe79PYn/fjCMehEhsGgRrFrl/G6J3RjjJaWWZUQkGshV1QwRqQb0AF4s0qYT8A4Qr6p7yiXSsnD1mvMzMng1ohXjlhzi4pj6vD2wB5E1q0DjxpCZ6fxejLOjavDl/V154/sNjJu3kYUb0xm0cjbXJKXQqFcv+OabkpOxq3xz7J8jmZjTgFfaDqCehPPxxi/o9tV7zjZWhjHGlJNSyzIicj7wIRCKs6c/SVVHicgoIEVVp4nI98B5wC7X235X1d4l7bdcyzLx8WTP/ZGhd4xhZnQb+p9VhVGDryAizPWHSuGSiQe95eTN6bw4ay2//p4BwLlpWxlwdBt9Jz5PjSrFfz4u2bafZ6auZM2uA/Rs15AX+pxPveVLYOhQZwMrwxhjysjTskyZau7eVC7J3ZW0M66/icGr4JeaZ/LMtW0YdOnZp3dzadd+1w19mgV7cpm+eAvLapxBraph9L+oKdeefybnNa5DaIiQdvAo363ezayVu/lpw14a1alKQq+29Gx3ht3g2hhz2jxN7sG1cFhiIjuSfuXO1rfxe+36vLV5Fr3C6oOcc9r7ZfZszgXOnTWLwbfDr7/vZ+LPW3h/4Vbe/WkLtaqEUZU89mbnoyLERFbn4StbMqTqXmo89jeP/0owxhhvCKrkvuaRp7mr7U4OV6/Fh2u/Iu6r9yFjw+nXtQsPf3T14i9ISOCCAXHsO5TDTxvS+GXRGvJmzuGMvTvoGR3CuV9/4uypxz/kvLCaklJ6nd4YY7wkaJL72t0H6J98hGpRUXx598Wcu6UeZO3wzlT+wuPZu3SB5GTn0gWLFlG/RgTXd2zM9U8NcibxyEh45Rvn6JfEROjTx5nY09PtAqoxpsIERXLfuvcQt09cTLXwUCbfF0fT+tXhDA8nGHlL0clN8fHHlyb45pvjF3CNMaYCBHxy352ZzcCJyRzLy+fTe12JvTyNHXtionY43I9+OZWZrMYY4yUBPVpm/6Ec+r3jYGfGET4d3IUOTeuWeZjjaSvcQ+/Z05K4MaZcBf1omayjx7jrg1/Ytu8wH/7tYmdih9NeqbHMEhKc9feCx8YY4wcCMrln/5zEkM+Ws7J2U8bfHktc88jjGwuXQypCXJzz4qkxxviRgEvuqsqj7y8kqUFbxq6eylVtrz2xgdW3jTEm8JI7wAUHttNxpYMbQ3b6OhRjjPFLAXcPVRFh0KO3Mrj+YefoFGOMMScJyJ67lV6MMaZkAddzN8YYUzpL7sYYE4QsuRtjTBAqNbmLSFURWSwiy0RklYgkumlTRUS+EJGNIpIsIjHlEawxxhjPeNJzPwpcoaodgI5AvIh0KdJmELBfVVsAYylyGz5jjDEVq9Tkrk5Zrqfhrq+iC9Jcj/NWfABfAleK3XbIGGN8xqOau4iEishSYA8wR1WTizRpDGwHUNVjQCYQiTHGGJ/wKLmrap6qdgSaABeLSPsiTdz10k9ablJEhohIioikpKWllT1aY4wxHinTJCZVzRCReUA8sLLQplSgKZAqImFAHWCfm/dPACYAiEiaiGzz8NBRwN6yxBok7Lwrl8p43pXxnOH0zvssTxqVmtxFJBrIdSX2akAPTr5gOg24E3AAfYEftJSF4lU12pMAXTGkeLJ+cbCx865cKuN5V8Zzhoo5b0967o2AD0UkFGcZZ5KqTheRUUCKqk4DJgIfi8hGnD32/uUWsTHGmFKVmtxVdTnQyc3rIws9zgZu9m5oxhhjTlWgzFCd4OsAfMTOu3KpjOddGc8ZKuC8fXYPVWOMMeUnUHruxhhjysCvkruIxIvIOtcaNU+52R6Ua9h4cN7DRGS1iCwXkbki4tFQKH9X2nkXatdXRFREAn5UhSfnLCL9XD/vVSLyaUXHWB48+B1vJiI/ishvrt/zv/oiTm8SkfdFZI+IrCxmu4jIm65/k+UicoFXA1BVv/gCQoFNwDlABLAMaFukzQPAeNfj/sAXvo67gs77cqC66/H9leW8Xe1qAQuARUCsr+OugJ91S+A3oJ7reQNfx11B5z0BuN/1uC2w1ddxe+G8/wJcAKwsZvtfgZk4J4F2AZK9eXx/6rlfDGxU1c2qmgN8jnPNmsKCcQ2bUs9bVX9U1cOup4twzhQOdJ78vAFGAy8B2RUZXDnx5JwHA+NUdT+Aqu6p4BjLgyfnrUBt1+M6QMDfIFlVF+BmMmch1wMfqdMioK6INPLW8f0puf+5Po1Lqus1t200eNaw8eS8CxuE89M+0JV63iLSCWiqqtMrMrBy5MnPuhXQSkQWisgiEYmvsOjKjyfn/SwwUERSgRnAgxUTmk+V9f9+mfjTPVQ9WZ/GozVsAozH5yQiA4FY4LJyjahilHjeIhKCc/nouyoqoArgyc86DGdppjvOv9B+EpH2qppRzrGVJ0/O+1bgA1V9VUTicE6KbK+q+eUfns+Uaz7zp557wfo0BZpw8p9mf7YpaQ2bAOPJeSMiPYCngd6qerSCYitPpZ13LaA9ME9EtuKsSU4L8Iuqnv6O/09Vc1V1C7AOZ7IPZJ6c9yBgEoCqOoCqONdfCWYe/d8/Vf6U3H8BWorI2SISgfOC6bQibQrWsAEP17AJAKWet6s88Q7OxB4MNVgo5bxVNVNVo1Q1RlVjcF5r6K2qKb4J1ys8+R2fivMCOiIShbNMs7lCo/Q+T877d+BKABFpgzO5B/vSsdOAO1yjZroAmaq6y2t79/UVZTdXj9fjvLL+tOu1UTj/U4PzBz4Z2AgsBs7xdcwVdN7fA38AS11f03wdc0Wcd5G28wjw0TIe/qwFeA1YDawA+vs65go677bAQpwjaZYCV/s6Zi+c82fALiAXZy99EHAfcF+hn/U417/JCm//ftsMVWOMCUL+VJYxxhjjJZbcjTEmCFlyN8aYIGTJ3RhjgpAld2OMCUKW3I0xJghZcjfGmCBkyd0YY4LQ/wO/wubGGOsS0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xArr, yArr = loadDataSet('ex0.txt')\n",
    "yHat = lwlrTest(xArr, xArr, yArr,0.01)\n",
    "xMat = mat(xArr)\n",
    "srtInd = xMat[:,1].argsort(0)   # 对输入变量按照第2列升序排列\n",
    "xSort = xMat[srtInd][:,0,:]\n",
    "\n",
    "\"\"\" 绘制图形 \"\"\"\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(xSort[:,1], yHat[srtInd])   # 拟合曲线\n",
    "ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0], s=2, c='red')   # 原始数据点\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "岭回归算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" 岭回归算法 \"\"\"\n",
    "def ridgeRegres(xMat,yMat,lam=0.2):\n",
    "    \"\"\" \n",
    "    输入：1、xMat——输入变量\n",
    "                2、yMat——目标值\n",
    "                3、lam——系数\n",
    "    输出：ws——回归系数\n",
    "    \"\"\"\n",
    "    xTx = xMat.T*xMat\n",
    "    denom = xTx + eye(shape(xMat)[1])*lam\n",
    "    if linalg.det(denom) == 0.0:\n",
    "        print(\"This matrix is singular, cannot do inverse\")\n",
    "        return\n",
    "    ws = denom.I * (xMat.T*yMat)\n",
    "    return ws\n",
    "\n",
    "\"\"\" 在1组lam上进行岭回归测试 \"\"\"\n",
    "def ridgeTest(xArr,yArr):\n",
    "    \"\"\" \n",
    "    输入：1、xArr——输入变量\n",
    "                2、yArr——目标值\n",
    "    输出：wMat——回归系数\n",
    "    \"\"\"\n",
    "    xMat = mat(xArr); yMat=mat(yArr).T\n",
    "    yMean = mean(yMat,0)\n",
    "    yMat = yMat - yMean     #to eliminate X0 take mean off of Y\n",
    "    # 正则化X（减均值，再除以方差）\n",
    "    xMeans = mean(xMat,0)   # 按列计算均值\n",
    "    xVar = var(xMat,0)              # 按列计算方差\n",
    "    xMat = (xMat - xMeans)/xVar\n",
    "    numTestPts = 30                # 测试次数\n",
    "    wMat = zeros((numTestPts,shape(xMat)[1]))    # 回归系数矩阵（1行代表1组系数）\n",
    "    for i in range(numTestPts):     # 每次变更lambda\n",
    "        ws = ridgeRegres(xMat,yMat,exp(i-10))\n",
    "        wMat[i,:]=ws.T\n",
    "    return wMat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在鲍鱼数据集上测试岭回归算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8XHd97//X55xZNNq30eLd8pLYsR3bUZyQOGQpoUlY0pQQQn+UUkpzaaEtl9vfD/q793eh/NpbbvtruW2h0GyFQFvgFggBEsjmBQhJLMdOYseJLduyLVu2JMvapdnO5/fHjMayLK8aaTQznyecx1nnnO/Jkd/zne985xxRVYwxxhQWJ9sFMMYYM/Ms/I0xpgBZ+BtjTAGy8DfGmAJk4W+MMQXIwt8YYwqQhb8xxhQgC39jjClAFv7GGFOAfNkuwLnU1tbqokWLsl0MY4zJKdu3b+9W1fCFtpu14b9o0SJaWlqyXQxjjMkpInLoYrazZh9jjClAFv7GGFOALPyNMaYAWfgbY0wBsvA3xpgCZOFvjDEFyMLfGGMK0Kzt53+5YpEEr/zsorq5Ti+Z5l3KmQeQSY6XXCbpF4/fRlIzIoI4qXlJzQuIkxqnlruu4LgOjivJwZecPnO5g8/v4C9yCYZ8uH4nfRxjzOySd+EfjyZoeaotu4WwxyID4DiCP+QSKPIlh5BLIDQ27aO8pojK+mIq64qpCIdw/fZB1JiZknfhHyoL8Imv3pbtYkwr1QnvLjrJpOq46dMrFD09ral9pcZj8+qdudzzFPWURFzxEoqX8NLjRGJsmeLFPeIxj9honMhInOhoglhqHB2NEx2NM9wXpff4MJHhOKNDsXS5RaCspij5RpB6Q6isD1FZV0xpdRGOY58gjMmkvAv/QnBWU4pMNjn7wzIyEqevc5jeE6mhc4TeE8N07O8gFkmktwuV+bni+kZW3thIVUNJFktsTP6w8DdZEwz5qFtYTt3C8jOWqyrD/VH6Ooc5dXyYw7t7eO25I+x85jANTRWs3NjIkvV1BIrsz9eYyyVnNSHMEs3NzWo3djNjhvujvPliB3t+2UHviWH8QZdlzXWsuHEO9YvL7YtlY1JEZLuqNl9wOwt/k0tUleP7+3jjhQ5aW04Qj3pUzylhxQ2NXHF9A6HSQLaLaExWWfibvBcdibOv5QR7XujgxMF+HJ/w9g8s56qb5ma7aMZkzcWGvzWampwVCPm46qa5XHXTXE4eHeSF77ey+V/fIjqSYN07F2S7eMbMatax2uSFmrml3PUHa1h6TR0vfL+Vl544cHaXWGNMmtX8Td5wfQ63/95V+IMuLU+2ERtNcOP7l9qXwcZMYso1fxGZLyKbRGSPiOwWkT+ZZBsRkX8QkVYReU1E1k/1uMZMxnGEWz90JWtuncerzx9h87fexPPsE4AxE2Wi5h8H/ouqviIiZcB2EXlGVd8Yt82dwLLUcB3w1dTYmIwTR9h43zICIR8tT7YRjSR4x++uxHWtldOYMVMOf1XtADpS0wMisgeYC4wP/7uBxzTZCPuiiFSKSGPqtcZknIhw3Xub8AddfvWD/cQjCX79gVX4/G62i2bMrJDRqpCILALWAS9NWDUXODJuvj21zJhptf7XF3LzB5fTtuskP/7ya0RH49kukjGzQsbCX0RKge8Bn1LV/omrJ3nJWQ2xIvKAiLSISEtXV1emimYK3Kqb5/GOj6zk2L5envj7nWfcUM6YQpWR8BcRP8ng/1dV/f4km7QD88fNzwOOTdxIVR9U1WZVbQ6Hw5komjEAXHFdA3f8/iq6Dg/w+Jd2MNwfzXaRjMmqTPT2EeARYI+q/t05NnsC+HCq18/1QJ+195uZ1rQuzLv+cA19J4b58ZdfRa0XkClgmaj53wj8NnCbiOxMDXeJyMdF5OOpbZ4EDgCtwEPAH2bguMZcsgVX1XDLh66k6/AAra90Zrs4xmRNJnr7/IIL3Dw+1cvnE1M9ljGZsOzaerY/1ca2n7SxZH2dPSjGFCTr+GwKjuMI1757Mac6hti/3Wr/pjBZ+JuCtHR9HdVzStj2k4P2C2BTkCz8TUESR7j2XYs5dXyY1pYT2S6OMTPOwt8UrCXrwtTMLWHbT9rwEl62i2PMjLLwNwVrrPbfe2KYfdus9m8Ki4W/KWhNa8PUzCu12r8pOBb+pqCJI2x492L6ukbY+7LV/k3hsPA3BW/x1bXUzi9l25NW+zeFw8LfFDyRZO2/v2uEt146nu3iGDMjLPyNARatqSW8oIyWJ9tIWO3fFAALf2MYV/vvHuWtX1nt3+Q/C39jUhaurqFuYar2H7fav8lvFv7GpIgIG97TxEDPKG/+yu44bvKbhb8x4yy4qpr6xeW0PGW1f5PfLPyNGWes7X+wJ8KeF6z2b/KXhb8xE8xfWU1DUznbn2ojEbPav8lPFv7GTJCs/TcxeCrCG78861HTxuQFC39jJjFvRRWNSyvY/lQb8Vgi28UxJuMs/I2ZhEjyaV9DfVH2bbOnfZn8Y+FvzDnMu6KKsuoi9u+w8Df5x8LfmHMQEZrWhznyRg+RkXi2i2NMRln4G3MeS9bV4SWUtte6s10UYzIqI+EvIo+KSKeI7DrH+ltEpE9EdqaG/56J4xoz3RoWl1NSEeDAjq5sF8WYjMpUzf/rwB0X2Obnqro2NXwhQ8c1ZlqJIzStDXN490liEev1Y/JHRsJfVbcCPZnYlzGzzZL1dcRjHod2ncx2UYzJmJls83+biLwqIk+JyFUzeFxjpqRxWSWhMr/1+jF5ZabC/xVgoapeDfwj8PhkG4nIAyLSIiItXV3WxmpmB8cRFl8d5tDrJ+0HXyZvzEj4q2q/qg6mpp8E/CJSO8l2D6pqs6o2h8PhmSiaMRdlybowsUiCI29Y66bJDzMS/iLSICKSmt6QOq41oJqcMfeKKoLFPvZbrx+TJ3yZ2ImI/DtwC1ArIu3A5wA/gKp+DbgX+AMRiQMjwP2qqpk4tjEzwfU5LF5Ty4FXu0nEPVyf/UTG5LaMhL+qfvAC678MfDkTxzImW5rW1/Hmi8dpf+sUC6+qyXZxjJkSq74Yc5Hmr6jCH3Q58Ir1+jG5z8LfmIvk87ssWl3DgVe78RL2kBeT2yz8jbkETevqGB2Mcay1L9tFMWZKLPyNuQQLV9Xg8zvst6Yfk+Ms/I25BP6gy4JVNRzY2YV61mHN5C4Lf2Mu0ZJ1YYb7ohw/YE0/JndZ+BtziRatrsXxif3gy+Q0C39jLlEg5GPBimr27+jEfqtocpWFvzGXoWldHYM9EToPDWS7KMZcFgt/Yy7D4qtrcRzhgN3m2eQoC39jLkNRiZ+5V1Sy/5Uua/oxOcnC35jL1LSujr6uEU4eHcp2UYy5ZBb+xlymprVhRLAffJmcZOFvzGUqLg/QuLTSunyanGThb8wULFkf5lTHEKeOW9OPyS0W/sZMQdPaOgD2v2K1f5NbLPyNmYLSqiD1i8vZb10+TY6x8Ddmipasr6P7yCB9XSPZLooxF83C35gpWrIuDGC1f5NTLPyNmaLy2hDhBWUcsF4/JodY+BuTAU3rwpw42M/gqdFsF8WYi2Lhb0wGnG76sdq/yQ0ZCX8ReVREOkVk1znWi4j8g4i0ishrIrI+E8c1ZraoaiihqrHEmn5MzshUzf/rwB3nWX8nsCw1PAB8NUPHNWbWWLIuTEdrL8P90WwXxZgLykj4q+pWoOc8m9wNPKZJLwKVItKYiWMbM1ssWR9GFQ6+arV/M/vNVJv/XODIuPn21LIziMgDItIiIi1dXfYPyOSWmrmllIdD1u5vcsJMhb9Msuysm6Cr6oOq2qyqzeFweAaKZUzmiAhL1oU5+uYpRodi2S6OMec1U+HfDswfNz8PODZDxzZmxjStC+N5Stvr3dkuijHnNVPh/wTw4VSvn+uBPlXtmKFjGzNj6heWU1oVtBu9mVnPl4mdiMi/A7cAtSLSDnwO8AOo6teAJ4G7gFZgGPjdTBzXmNlGHKFpbZjdPz9GdDROoCgj/8SMybiM/GWq6gcvsF6BT2TiWMbMdkvWh3ltUzuHdp1kWXN9totjzKTsF77GZFjDkkpCZX77wZeZ1Sz8jckwxxEWrw3Ttusk8Wgi28UxZlIW/sZMgyXrwsQjCQ6/cb7fPhqTPRb+xkyDuVdUESz2cWCnNf2Y2cnC35hp4LoOi9bU0vZaN4m4l+3iGHMWC39jpsmSdWEiw3GOvnUq20Ux5iwW/sZMk/krq/EHXfZb04+ZhSz8jZkmPr/LwtU1HNzZheeddSsrY7LKwt+YadS0NszIQIyO1t5sF8WYM1j4GzONFq6qwfU7dptnM+tY+BszjQJFPhasrObAji7Umn7MLGLhb8w0W7IuzFBvhBOH+rNdFGPSLPyNmWYLV9fiOMIBu82zmUUs/I2ZZkUlfuZdWcX+HZ0kb3BrTPZZ+BszA5rWhenvHuXk0cFsF8UYwMLfmBmx+OowItgTvsysYeFvzAwoLg/QuLTSunyaWcPC35gZsmR9mFMdQ5w6PpTtohhj4W/MTGlaWwdgtX8zK1j4GzNDSquC1C8ut8c7mlnBwt+YGdS0LkzX4QH6u0eyXRRT4DIS/iJyh4i8JSKtIvLZSdZ/RES6RGRnavhYJo5rTK5Zur4OBHZtPZrtopgCN+XwFxEX+ApwJ7AS+KCIrJxk0++o6trU8PBUj2tMLiqvDbGsuZ7XN7czMhDNdnFMActEzX8D0KqqB1Q1CnwbuDsD+zUmL137rkUkYh47njmc7aKYApaJ8J8LHBk3355aNtH7ROQ1EfkPEZmfgeMak5OqGkpYdm2y9j/cb7V/kx2ZCH+ZZNnEG5j8CFikqmuAZ4FvTLojkQdEpEVEWrq6rEeEyV/NdyVr/zut9m+yJBPh3w6Mr8nPA46N30BVT6pqJDX7EHDNZDtS1QdVtVlVm8PhcAaKZszslK79b7Hav8mOTIT/NmCZiCwWkQBwP/DE+A1EpHHc7HuBPRk4rjE5baz2b23/JhumHP6qGgc+CfyMZKh/V1V3i8gXROS9qc3+WER2i8irwB8DH5nqcY3JdVUNJSzbUM8ua/s3WSCz9f7izc3N2tLSku1iGDOtek8M82+ff5Gr37GAG9+3NNvFMXlARLaravOFtrNf+BqTRZX1xSzf0GC1fzPjLPyNybLmuxaRiHvsePpQtotiCoiFvzFZlq79bzlqtX8zY3zZLoAxJln73/vycbY91cryt1fR09NDT08P0WiUYDBIUVERwWBw0ulAIIDjWD3OXBoLf2Nm2PDwMCdPnuTUqVPpkO/p6eFUYxeb9kTYNK4jtIhc1EPfi4qKmD9/PkuXLmXp0qXU1NRM4xmYfGDhb8w0SyQStLe309rayr59+zh+/PgZ68vLy6muruaKK67g4MsDNK2cxw13raS6uppAIEAsFmN0dJRIJEIkEpl0enBwkIMHD7Jv3z4AqqqqWLZsGUuXLmXRokUEAoFsnLqZxSz8jZkGAwMD6bA/cOAAo6OjiAjz58/ntttuo76+nurqaiorK/H7/enXPTvyBvu3d1J+bzXBYBCAQCBw0eHd09OTPu4rr7zCyy+/jOu6LFy4kKVLl7Js2TJqa2sRmeyuLKaQWD9/YzJgrHa/b98+Wltb07X70tLSdA28qamJUCh03v30dg7zb59/iTW3zWPjvcumVKZYLMbhw4fTbwbd3d0AhMNhNm7cyKpVq3Bdd0rHMLPPxfbzt/A35jINDw/T2trK3r17aW1tPaN2Pxb4DQ0Nl1zLfu7rb9C6vZMP/cXbKKkIZqy8vb297Nu3j23bttHZ2UlVVRU33XQTa9asweezRoB8YeFvTIapKt3d3ezdu5e9e/dy+PBhVJXi4mKWL1/OsmXLLqp2fyHp2v+t89j4/qnV/ifjeR5vvfUWW7dupaOjg4qKCjZu3MjatWvPaIIyucnC35gMiMfjHDp0KB34p06dAqC+vp7ly5ezfPly5s6dm/Guls994w32tXTy2xmu/Y+nqrS2trJlyxba29spKyvjhhtu4JprrrEviHOYhb8xlyiRSNDV1UVHR0d6OH78OLFYDNd1aWpqSgd+RUXFtJZlrPa/6u1zefv9y6f1WKrKwYMH2bp1K21tbRQXF3PDDTdw7bXXpr90NrnjYsPfGvpMQYrFYnR2dp4R9CdOnCCRSADg9/tpaGhg3bp1LFmyhMWLF89obbiyrpiVG+fw+uZ2ikr9XPuuRdPWQ0dEaGpqoqmpiUOHDrF161aeffZZfvnLX/K2t72N6667zt4E8pDV/E1eicViDA4OpoehoaEz5seGvr6+9I+nioqKaGxspLGxkYaGBhobG6mpqcn6r2YTCY/N33yTN188zqqb53LTB5bjODPTRfPo0aNs2bKFvXv3EgqFuOGGG9iwYYO9CeSAgm32icfj7N+/fxpKlH8u9tpP3G78/GTrxpZdaFpV8TzvnMPY+ng8TjweJxaLnXOIx+NEo1Gi0cnvjVNcXExpaWl6qKioSAd+ZWXlrO33rqr86gf72fH0YZasD3P7716F65+5N6WjR4+yefNm9u3bl24O2rBhg30nMIsVbPgPDQ3xN3/zN9NQIjPTHMdJD4FAAL/fj8/nw+/3n3MoKSmhpKTkjKAvKSnJ+f7sO545zAvfa2XuFVXc9fHVBEIz22Lb3t7Opk2b2L9/P8XFxWzcuJHm5mZ7E5iFCjb8R2OjvLz35WkoUe6Yci12kpfL2EI5/7L08SW1fvw4Va6x+9WIIziOg4ggjiCSmk9Nj1HO/huddNnETyHoWcvHL5s4rejp5Up62cTtx+9z7DVjHEmejyMODuOmxUmvd8Qh4AYIOkGCviBFbhEBN0CRr4iAEzjn9XvrxQ6ef+xNauaV8u5PXk1x+cwH75EjR9i0aRMHDhygpKQk/SZgXURnj4IN/57RHm7+zs3TUCJjZkbQDSbfDNwiiv3FVAYr00N15wKCzzfhligLPiiEGyqoClZRWVRJdVF1+k1muh06dIjNmzdz8OBBSktLuf7661m1ahWVlZUzcnxzbgUb/tFElBc7XpyGEplskkk+jkxWQx7bbuKnkvGvT38CGfvf+HmRM8aTbT+2bPwxBEl/YvDUwyP1nYV6eOqhaHqc8BJEE1EiXoRIPEIkMcmQWj4UG6I30ktfpI9TkVP0RfooP1XPnXseIOHEeXLF1zhZcgwAn+OjvriexpJGGksaaShpoKGkIT3fWNpIib9kKpfhLG1tbWzevJm2tjYA5s2bx1VXXcXKlSunvTusmVzBhr8x+W40PsrhQyfY+mAbsdEEDe+LMVrXw4mhE3QMdXB86DgdQx10DneS0MQZry0LlNFY0sic0jnMK53HnNI5Z0yXBcouq0w9PT3s3r2b3bt3p+9rtGDBgvQbQVnZ5e3XXDoLf2Py3EDPKD/6h530d4/yzt+7iqZ14TPWx7043SPd6TeDjqEOjg0e4/jQcY4OHuXo4FFG4iNnvKYsUHbGm0JDcQP1JfXUF9fTUNJAbagWn3P+L5u7u7vTbwSdnZ0ALFy4kFWrVrFixQpKS0sz+x/CnMHC35gCMDoY48dfeZXOtn7mrahm0eoaFq2upbz2wvcXUlV6I70cGzzG0cGj6fHY9LGhY2e9OTjiUFtUm35DqC+pp6G4gdriWqqLqqkpqqG6qJrKokr8jp/Ozs70G8HYXUWrq6upq6ujrq6O+vp66urqqK6uzvkeWbPFjIa/iNwB/D3gAg+r6hcnrA8CjwHXACeBD6hq2/n2ebnh3zcS4+Pf3H7JrzMzb8qdks7TK2niuvHfD0h6GTiSfEVytZDqqJQaJ+cdR3BEcFPbJ+fBTS13RNLTfp/gdxx8ruB3Hfyu4HMc/D4HvyP4UstCfpfigI9QwKU4NSSnfYT8Lu4l/JgrFknQ8mQbB3Z20XtiGIDqOSUsWlPL4jW11C0qv6wfh6kqfaO9nOhtp6vvGF39HZzsP05PfyenBjvpHeiib7CbRGQUX4Lk4JGeLpMQ5U4x5RKiREL4nXJGpZJRp4hhCTCMm75QDkolUINQDVR7SoUIxT5f8jGVgQDi9yeHidMBP05JCW55OU5pKW55OW5ZGVJcPGt/vzGdZuz2DiLiAl8BbgfagW0i8oSqvjFus98DTqnqUhG5H/ifwAemeuzJaCTC6td/Ph27Nhk0WVfNDOz0jH3rhOXjJ0WTW2mqS2f6tTrxtYqniurY9Oll8XHTCqiX/FFawjvd/TP9RoMi4wokqTKcXqepbZNjn0DAFYp8QrHPIeRzCPmEkE8och1CLgR9Qsh1CLhQ4whzHCVSE6RjtIrjnTXs+OkAr/z0EAEi1Mtx6r12woljuPFRNBZLDtHo6emJQzQK8TgANanh0gylhjMlBDwH4q5LX3kZfZWV9FdUMFBeweHyCvaFimHszSoRxxmMEBwdpSgSoWh0lODYeDSSnI6M4ovF8cdi+OJxfPHktCOCO/amUFaGW1aGU16GW1GBW16RHFeU41ZU4Iwtq6xIby95/lzkTPxSZAPQqqoHAETk28DdwPjwvxv4fGr6P4Avi4joNLQ5lXpR7nnm0Uzv1phZwUNQETwET5LTMXGIIiTEIeE4+ByXua5LfaCEgcor6K9YzrHSJRxxFyJunBKvG18ggs8fwV8cwa8R/Izg1xECOkxABwnoIP5EP0IMcYDUkPxNxth8alpOzyPJT0vqJMekyiiOD5EylHI8pwiVIEqAEgkSIkCDBIEAQoB4zGUIjxHxiBAnInEiRXGioQTDxOmVGBHieHL++BAFnzq4KvjGBk9wEuCejOF2deN4XTiquJ7iJDwcTY09D0m+pad6c6XGounzFscB10lNC+K6qf8WDuI66fWO44Dr4rgOOA7qpD5Vpr47EdeB1CfH5H9Dh1BFOTfc9/5p+zuCzIT/XODIuPl24LpzbaOqcRHpI1mR6M7A8c/gVlSw5NlnM71bkyMu6VP+pO1GF1gmp+vyp5eNXy2ntznXePx24wdOB+b4QZxkaJCIINEhEpFBhgZOMdjfy/BgH6OD/USG+5ChHnzDPcjoKfyRXoKxXhbFt1Ca+BGl3jAnY0s5GLmWU/H5xDREVEMMepVENURUS1Auv6ZbpBBSISQQktPjYgeKHKFoYrOTnh6rKjH1iGucuBcjplGKvFESXhTFB6kusuCNe5lHHI+oeEQdj4SjyU8UAgkHEpKcT4gSFyUhXnLsJscRFE+UBB4JUTw84ggqkGy9nsr3D15yUCB+eXuoOhTihvumUISLkInwn+yf28S35IvZBhF5AHgAkt3ELqswsUECWz51Wa81M20aGv2ZGLRnB++Zr5dLHDNhGWeu01RtUUmNz0i5ZN//aBGJqIvGXLyYixf34cV8aNyPFw/gJQJ4iSCeVwSeIDqK6AhCBCGKSBSIUkyUEiKIxBAiODKI4wzjBAWn2IdTXodTthAproTiauqDVZT5yhl2Kxj1VzDsq2TIKWfYLWM4LoyOxhkdjhMZiTE6Eic2mkAV3IRSHPMIxTxCUY/iWGqIJodQ7Oy3jZgDI0GXkYBDb9Bh2OcxMtzJUP9hIiMniQ33EBvqJhYfJqGx0y/0BXBKqnBKKnFKyxBfAFw/4voRnx9xA+BLzY8tc3zgnA5rRVIV9vHXXPAp+FSTHwm8sXY+LzmtiniKxOOIF0MSMVAP8RKgcdTzUPWQ1Bj1Upc0+QlB1Us266XGSqpZb9zlH9cumf7LVFKfKias9vm8s/9eMywT4d8OzB83Pw84do5t2kXEB1QAPRN3pKoPAg9C8gvfyyqNejDUdVkvNTNoyi1+k7xe9cx1kzb6nxnGlzRO7+pc6zT9ppDQUuKJeuKJhuQ4Xkc8UUc8Xody7tsyiBPBcaM4bhQJxBAXVCvx1I+qLzl4Dppw0GTV9uydRID+cfsMusk3g2I/TsiHz+9QqlDiKWGNgteNepqssKompxU07uENRPGGJ1RfHcGtCODWhfBVFuFWBtODryI5lqDLyfbDHNzRwtEd2zj21h68RIJAqJjqOXOpnd9AWe1qysN1lNWGKa+to7w2TFFpWUF+SZsNmQj/bcAyEVkMHAXuB35rwjZPAL8D/Aq4F3h+Otr7AQhVwQObp2XXxoynCY/4yVHincPEuoaJd40Q704OZwSmI/iqi/DVhgjWhvDVFuGUBHBCLk6RDyeUHCToQ9xLCz5VhbiisQTeSDw5DMfxhmNnTg+PrYvhDcWSbctj7fJjbfU+wHGSPYNS7c/u4opkqFcGcatSQV8WQCbpPRSLjHJk9+sc+PE2Du5sob8r2cc/vGARze/5TRava2bOsitxrEvnrDDl8E+14X8S+BnJhrJHVXW3iHwBaFHVJ4BHgG+KSCvJGv/9Uz2uMTPFiySIdw0T6xpJBn3nMPHOYeInR083HwBueQBfOERodS2+2mJ84RC+2hC+qmDyS71pICLgF8Tv4BTP/M3VVJW9L/6C3Zuf5cju14nHoviCQRauXst1v3Efi9c1U1ZTO+PlMhdmP/IyBc2LJEj0R0j0R5NDXwSvP3p6WV+ERN+4ZwQ44KsJ4QsX468rxlcXwh9Ojp1gYT0Y79jeN9n82EN07HuLivoGlqzfwOJ1zcxbuRqf3eUza+wxjibn6NgXZKc74CeXeYomxsYeJMbNx71kG3UiuU4jHl40gUbiyelIHI0k8CIJNJpIT3vDMRJ9UTSSOKscEnRxywO4FUGCSyrx1YaSQR8O4asJIb787v99If1dnWz9t6/z1gtbKams4tc//iesvPk2HMeac3JJ3oV/YijGiS8V8C98L+mD3EVufJ7NzvrgmP5C9PzLdPyy8aGfaQIScJNfegZT44CLP1xM0dIqnPJAMujLg8kvMcsDBVeDv1jRkWFeevx/s/0njyMI17/vfq597/sIFF34VhJm9sm7v3LxCaGrLv23iHllOnpLnGeXk/bOOKOX3Zm3Vx7rHnm6TzvpeyqklzkT5l1J/nDGTX1R6Z5ehpt6GIwr44LehwRdxO9M+uWkuXiel2DXpmf45Xe+xXBfLytuupWN93+Y8trwhV9sZq28C38n6KPqnmXZLoYxeeHQazvZ/M2H6T7cxpwrVvIb/9f/Q+PSK7JdLJMBeRf+xpipG+o9xdMP/iMHtr9Mebied3/qsyy//kZbC08hAAANR0lEQVTrg59HLPyNMWfoaH2LJ/72fzA6OMhNv/UR1t/5Xnz2oPa8Y+FvjEnbtekZnn3knyiprOaD/+/fULeoKdtFMtPEwt8YQyIeZ/NjD7PzZz9mwaqrefenPkOorDzbxTLTyMLfmAI33NfLj770Rdr37OKad9/D23/rI3YLhgJg4W9MATu+fx8//Nu/ZHRggLv+6E9ZsfGWbBfJzBALf2MK1O4tz/HMQ1+mpLKK+7/w19QvXpLtIpkZZOFvTIFJxONs+dYj7HjqR8y/ag3v/tRnKC6vyHaxzAyz8DemgAz39/GjL/0V7W/sYv1dd3Pzhz5q7fsFysLfmAJxvHUvT/zdXzHS38edn/g0K99+W7aLZLLIwt+YAvD680/z3CP/RElVdbJ9v2lptotksszC35g8Fo/F2PQv/8xrz/2UBavX8q4//j+tfd8AFv7G5K2Bk9088Xf/g+Ote9lw973ceP9v2z33TZqFvzF56Mgbr/Pj//U/iUUivOfTf8by627MdpHMLGPhb0weUVVeefKHbPnWo1Q2zOG+//5X1Mybn+1imVnIwt+YPBEbHeXpB/+RN3+5haXXXs8df/hpgsXF2S6WmaUs/I3JA73HO/jh3/4l3UcOsfH+D7Ph7nsRp7CfNWzOz8LfmBymquza9AxbvvkIIsL7Pvt5Fq29JtvFMjlgSuEvItXAd4BFQBtwn6qemmS7BPB6avawqr53Ksc1xsDJo0d49qGv0L5nF3OvXMmdn/g0FXUN2S6WyRFTrfl/FnhOVb8oIp9NzX9mku1GVHXtFI9ljCHZd//lx7/Ly4//b3zBILc/8EesvvV2a+Yxl2Sq4X83cEtq+hvAZiYPf2NMBhx543WeeegrnDrWzpU33swtH/4YJZVV2S6WyUFTDf96Ve0AUNUOEak7x3ZFItICxIEvqurjk20kIg8ADwAsWLBgikUzJn+MDPSz5VuPsnvzs1TU1fObf/bnLLa2fTMFFwx/EXkWmKwh8b9ewnEWqOoxEWkCnheR11V1/8SNVPVB4EGA5uZmvYT9G5OXVJU9v9jM5sceZnRwgGvvvpe3ve9+/MGibBfN5LgLhr+qvuNc60TkhIg0pmr9jUDnOfZxLDU+ICKbgXXAWeFvjDmt93gHzz7yTxx6bQeNS6/g9v/2F4QXLs52sUyemGqzzxPA7wBfTI1/OHEDEakChlU1IiK1wI3AX0/xuMbkrVhklJd/+B9se+J7uD4/v/bRP2DN7XfYfXlMRk01/L8IfFdEfg84DLwfQESagY+r6seAFcA/i4gHOCTb/N+Y4nGNyTuqSmvLi2z+xkP0d3Vy5Y03c/OHPkppdU22i2by0JTCX1VPAr82yfIW4GOp6ReA1VM5jjH57lTHUZ7/+oO07dxO7fyF3Pe5v2L+SvtnY6aP/cLXmCyKjY7y0uPfpeVH38f1B7jlw7/P2l9/F67P/mma6WV/YcZkgarS+vKv2PTYQwx0d7Hyplt5+4c+an32zYyx8DdmhvUcO8rz//I1Dr22g/CCRdz153/KvCuvynaxTIGx8DdmhnhegpYf/YAXvvstXH+AWz/yn1j7zrtwXOvFY2aehb8xM+DU8WP89Ctf4tjePSy77gZ+7aN/YE08Jqss/I2ZRqrKq888xZZvPYLr83HXH/0pV954MyKS7aKZAmfhb8w0GTjZzc++9vccem0Hi65ezzs//seUVddmu1jGABb+xmScqvLmLzbz3L98jUQ8zjs+9oesecedVts3s4qFvzEZNNzfx7MPf4V9L73AnOUruOMT/5mqhjnZLpYxZ7HwNyZD9m9/iaf/+R+JDA1y0299hOb33GP34zGzloW/MVM03N/H1m89yu4tzxFeuJh7/9tfEF6wKNvFMua8LPyNuUyqyhtbn2fzNx8hOjzEdfd8gLfdez+uz5/tohlzQRb+xlyGUx1Hefbhf+LwrleZs3wFt//+J6i12r7JIRb+xlyCRDzGtie+z4vf/zY+fyDZk+fX7rCHp5ucY+FvzEU6+uYbPPPQlznZfpjl12/k1o88QGlVdbaLZcxlsfA35gJGhwb5+b9+ndee+ylltWHu+cznaFp/bbaLZcyUWPgbcw6el2Dvr37B5sceZrivj2ve9RvccN//QaAolO2iGTNlFv7GTDDQ082uTc/w+vNPM9DdRX3TUu75zOeob1qa7aIZkzEW/sYAXiLBwZ0tvPbsTzm4YzuqHgtWr+XmD32UZRtusNsum7xj4W8KWn9XJ69veppdzz/N4Kkeiisqufbu97H61ndS2dCY7eIZM20s/E3BicdiHHxlG689/zPaXn0FgMVXr+e2j36cpvUb7Pm5piBM6a9cRN4PfB5YAWxQ1ZZzbHcH8PeACzysql+cynGNuVgjgwN0tR2k69ABOtsO0NV2gJNHj+AlEpRW13D9b97P6ltvpzxcl+2iGjOjplrF2QX8JvDP59pARFzgK8DtQDuwTUSeUNU3pnhsY4DkbRZioyMM9fXSfeQQXW0H6Gw7SGfbfga6u9LblVRVU7dwMU3XbGDOFStYtGa9teWbgjWl8FfVPcCF7lO+AWhV1QOpbb8N3A1Y+Bco9Tw8z8OLx4nHosRjURKxOIlYlHgsRiIWJRGLpaZjREdHGOnvZ3Swn5GBfkb6+xkZHEhOD/QzOtBPIh5P71/EoWrOXOZesZLwOxdTt6iJ8MLF9thEY8aZicbNucCRcfPtwHXTdbCRwQG+87nPTNfuL4qqZvPgZy+64LY6btG4rXVspKCa2lyT2+jp5agmN9XkOs/zIBXwmho8z0M1OX25RByKSksJlZUTKi+noq6BhiXLCZWXEyotI1ReQc28+dTOX4g/WHTZxzGmEFww/EXkWaBhklX/VVV/eBHHmOxjwaR5JCIPAA8ALFiw4CJ2fTbHcaiZO/+yXptR2Xxq0yTHPmdpUtue69NberlIch8iqWWS/H9qOjlKrhPHxXGc1LSDOE5yfmyQ5Lzr9+P6/fj8gfR0et7nxxfw4/r8+ItChMrLKSousXvoGJMhFwx/VX3HFI/RDoxP43nAsXMc60HgQYDm5ubLqj4Hi0t4z6f/7HJeaowxBWMmqlHbgGUislhEAsD9wBMzcFxjjDHnMKXwF5F7RKQdeBvwExH5WWr5HBF5EkBV48AngZ8Be4DvquruqRXbGGPMVEy1t88PgB9MsvwYcNe4+SeBJ6dyLGOMMZlj354ZY0wBsvA3xpgCZOFvjDEFyMLfGGMKkIW/McYUIMnqrQjOQ0S6gENT2EUt0J2h4swG+XY+kH/nlG/nA/l3Tvl2PnD2OS1U1fCFXjRrw3+qRKRFVZuzXY5Mybfzgfw7p3w7H8i/c8q384HLPydr9jHGmAJk4W+MMQUon8P/wWwXIMPy7Xwg/84p384H8u+c8u184DLPKW/b/I0xxpxbPtf8jTHGnEPehb+I3CEib4lIq4h8NtvlyQQRaROR10Vkp4i0ZLs8l0pEHhWRThHZNW5ZtYg8IyL7UuOcesbiOc7p8yJyNHWddorIXefbx2wiIvNFZJOI7BGR3SLyJ6nlOXmdznM+uXyNikTkZRF5NXVOf55avlhEXkpdo++kbp1/4f3lU7NP6mHxexn3sHjgg7n+sHgRaQOaVTUn+yeLyNuBQeAxVV2VWvbXQI+qfjH1Jl2lqtl9/uYlOMc5fR4YVNX/L5tluxwi0gg0quorIlIGbAd+A/gIOXidznM+95G710iAElUdFBE/8AvgT4BPA99X1W+LyNeAV1X1qxfaX77V/NMPi1fVKDD2sHiTRaq6FeiZsPhu4Bup6W+Q/IeZM85xTjlLVTtU9ZXU9ADJZ2/MJUev03nOJ2dp0mBq1p8aFLgN+I/U8ou+RvkW/pM9LD6nL3iKAk+LyPbUc47zQb2qdkDyHypQl+XyZMonReS1VLNQTjSRTCQii4B1wEvkwXWacD6Qw9dIRFwR2Ql0As8A+4He1EOz4BIyL9/C/6IfFp9jblTV9cCdwCdSTQ5m9vkqsARYC3QAf5vd4lw6ESkFvgd8SlX7s12eqZrkfHL6GqlqQlXXknwW+gZgxWSbXcy+8i38L/ph8bkk9WQ0VLWT5JPTNmS3RBlxItUuO9Y+25nl8kyZqp5I/eP0gIfIseuUakf+HvCvqvr91OKcvU6TnU+uX6MxqtoLbAauBypFZOypjBedefkW/nn3sHgRKUl9YYWIlADvBHad/1U54Qngd1LTvwP8MItlyYixkEy5hxy6TqkvEx8B9qjq341blZPX6Vznk+PXKCwilanpEPAOkt9lbALuTW120dcor3r7AKS6bv0vwAUeVdW/zHKRpkREmjj9nGQf8G+5dk4i8u/ALSTvPngC+BzwOPBdYAFwGHi/qubMF6jnOKdbSDYnKNAG/Kex9vLZTkQ2Aj8HXge81OL/m2Q7ec5dp/OczwfJ3Wu0huQXui7Jivt3VfULqYz4NlAN7AA+pKqRC+4v38LfGGPMheVbs48xxpiLYOFvjDEFyMLfGGMKkIW/McYUIAt/Y4wpQBb+xhhTgCz8jTGmAFn4G2NMAfr/AfEvICh2KzpGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xArr, yArr = loadDataSet('abalone.txt')\n",
    "ridgeWeights = ridgeTest(xArr, yArr)       # 岭回归算法计算出的系数\n",
    "\n",
    "\"\"\" 绘制图形 \"\"\"\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(ridgeWeights)   # 回归系数曲线\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "前向逐步线性回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" 按列正则化 \"\"\"\n",
    "def regularize(xMat):\n",
    "    \"\"\" \n",
    "    输入：xMat——输入变量\n",
    "    输出：inMat——正则化结果\n",
    "    \"\"\"\n",
    "    inMat = xMat.copy()\n",
    "    inMeans = mean(inMat,0)   # 按列计算均值\n",
    "    inVar = var(inMat,0)             # 按列计算方差\n",
    "    inMat = (inMat - inMeans)/inVar\n",
    "    return inMat\n",
    "\n",
    "\"\"\" 前向逐步线性回归 \"\"\"\n",
    "def stageWise(xArr,yArr,eps=0.01,numIt=100):\n",
    "    \"\"\" \n",
    "    输入：1、xArr——输入变量\n",
    "                2、yArr——目标值\n",
    "                3、eps——系数变化量（步长）\n",
    "                4、numIt——迭代次数\n",
    "    输出：ws——回归系数\n",
    "    \"\"\"\n",
    "    xMat = mat(xArr); yMat=mat(yArr).T\n",
    "    yMean = mean(yMat,0)\n",
    "    yMat = yMat - yMean\n",
    "    xMat = regularize(xMat)      # 正则化输入变量（满足0均值和单位方差）\n",
    "    m,n=shape(xMat)\n",
    "    ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy()   # 初始化所有系数为1\n",
    "    for i in range(numIt):\n",
    "#         print(ws.T)\n",
    "        lowestError = inf;         # 每轮迭代设置当前最小误差为正无穷\n",
    "        for j in range(n):          # 遍历每个特征\n",
    "            for sign in [-1,1]:     # 系数变化方式（增或减）\n",
    "                wsTest = ws.copy()\n",
    "                wsTest[j] += eps*sign                   # 微调系数\n",
    "                yTest = xMat*wsTest\n",
    "                rssE = rssError(yMat.A,yTest.A)   # 计算新的误差（×.A从mat转为array）\n",
    "                if rssE < lowestError:                    # 在新误差<当前最小误差时，保存系数\n",
    "                    lowestError = rssE\n",
    "                    wsMax = wsTest                        \n",
    "        ws = wsMax.copy()\n",
    "    return ws"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在鲍鱼数据集上测试前向逐步线性回归算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "前向逐步线性回归计算系数为： [[ 0.043 -0.011  0.12   0.022  2.023 -0.963 -0.105  0.187]]\n",
      "标准线性回归计算系数为： [[ 0.0430442  -0.02274163  0.13214087  0.02075182  2.22403814 -0.99895312\n",
      "  -0.11725427  0.16622915]]\n"
     ]
    }
   ],
   "source": [
    "stageWeights = stageWise(xArr, yArr, 0.001, 5000)   # 前向逐步线性回归算法计算出的系数\n",
    "xMat = mat(xArr); yMat=mat(yArr).T\n",
    "xMat = regularize(xMat) \n",
    "yM = mean(yMat,0)\n",
    "yMat = yMat - yM\n",
    "ws = standRegres(xMat,yMat.T)                                    # 标准线性回归算法计算出的系数\n",
    "print('前向逐步线性回归计算系数为：', stageWeights.T)\n",
    "print('标准线性回归计算系数为：', ws.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
